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Abstract 

We present a stochastic description of multiple scattering of polarized 
waves in the regime of forward scattering. In this regime, if the source 
is polarized, polarization survives along a few transport mean free paths, 
making it possible to measure an outgoing polarization distribution. We 
solve the direct problem using compound Poisson processes on the rota- 
tion group SO(3) and non-commutative harmonic analysis. The obtained 
solution generalizes previous works in multiple scattering theory and is 
used to design an algorithm solving the inverse problem of estimating the 
scattering properties of the medium from the observations. This tech- 
nique applies to thin disordered layers, spatially fluctuating media and 
multiple scattering systems and is based on the polarization but not on 
the signal amplitude. We suggest that it can be used as a non invasive 
testing method. 

1 Introduction 

Soon after its discovery in adiabatic quantum systems by Berry [3] , experimental 
evidence of a geometric phase for polarized light travelling along a bended optical 
fiber was shown by Tomita and Chiao [23 . It was rapidly realized that any 
classical transverse polarized wave could exhibit a geometric phase [52] when 
travelling along a curved path in a three dimensions space. This, however, is 
only possible if the wave remains polarized, and if the direction of propagation 
evolves smoothly with time. 

If one emits a plane monochromatic wave into a scattering medium, the 
multiple scattering events will result in a spreading of the distribution of the 
wave vector with time [25] . Depending on the strength and density of the 
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scatterers, this spreading will happen at different rates. The correlation length 
of the wave vector is called the transport mean free path £*. It is directly related, 
in the multiple scattering regime, to the effective diffusion constant of the energy 
in the system D — c£*/3 [TH] (c is the velocity of the wave). Several multiple 
scattering regimes exist, depending on the ratio between the wavelength A, the 
size a of the scatterers, the mean free path £, the transport mean free path £* 
and the depth of the medium L [9]. In this article, we will investigate the regime 

A < a « £ < t < L. (1) 

The condition £ <C £* states that the wave is preferentially scattered to a di- 
rection close to the incoming one. This is the so called "forward scattering" 
regime. For spherical scatterers, this is a consequence of the condition A < a. 

The transport mean free path £* is the length scale of disorder in the medium 
while the transport mean free path is the scale along which the waves behave 
as homogeneous. Measuring the transport mean free path is therefore an im- 
portant issue in experiments because it is the resolution scale for most of the 
physical quantities measurable using waves. Several estimation techniques al- 
ready exist for the transport mean free path that can be distinguished from the 
direction of the measured signal. The simplest idea is to use the intensity trans- 
mission ratio [7] to measure £* thanks to the diffusion approximation. Diffusion 
models are also useful in anisotropic media [TU] and to describe intensity leaks 
in direction orthogonal to the incoming source [14] . In backscattering setup, the 
opening angle of the weak localization cone, in which an intensity enhancement 
due to interferences is observed, is related to the transport mean free path in 
the medium |24j . All these experiments rely on intensity measurements. An 
amplitude independent method has been proposed to measure the mean free 
path £ from the phase statistics of a Gaussian field p] . 

It is known that in the forward scattering regime polarization is parallel 
transported during the propagation, which may result in the existence of a geo- 
metric phase |15) . In the eikonal approximation, each ray possesses a geometric 
phase depending on its geometry. At a given point, the observed phase is not 
uniquely defined, but has a distribution related to the distribution of paths from 
the source. This distribution was recently shown to depend on the outgoing di- 
rection of the ray [2U] . 

The distribution of geometric phase from multiple scattering of a polarized 
incident wave was adressed ten years ago [12] without any condition on the 
final direction. The calculation was based on Brownian motion at the surface 
of a sphere and followed an approach developed earlier by Antoine et al. [2]. 
However multiple scattering is rigorously a Brownian motion on a sphere only 
in the limit of very weak and dense scatterers. A more general model would be 
a random walk with macroscopic steps. 

We consider media with a thickness L of the order of a few £* , therefore the 
outgoing wave vectors are not evenly distributed. The joint distribution of di- 
rection and polarization state, related to paths statistics, contains informations 
concerning the scattering events that are not yet randomized. The number of 
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scattering events is of the order of Ljl and in the case where this number is not 
very large (L/l ~ 10) the usual approximations are not suited. For this reason, 
attention was recently brought to stochastic processes for the description of a 
scattering system where the number of scattering events is small. 

Stochastic models have been introduced fifteen years ago in multiple scatter- 
ing [IB] to solve direct problems. Recently Le Bihan and Margerin [T3] showed 
that a stochastic model for the direction of propagation can be used to solve 
inverse problems. They computed the mean free path in a medium from the 
distribution of outgoing directions. We extend this model to take polarization 
into account and estimate the transport mean free path I* from the geometric 
phase distribution. 

2 Scattering of polarized waves 

Without loss of generality, we consider linearly polarized waves. The results we 
present remain valid as long as polarization is not purely circular. 

Polarization of transverse waves lies in the plane orthogonal to the direction 
of propagation. We describe linearly polarized plane waves with two vectors: 
The direction of propagation k and the direction of polarization p, with p _L fe. 
We consider only media with isotropic uniform absorption, the wave amplitude 
is therefore completely determined for all paths and does not play any relevant 
role. The direction k lies on the unit sphere S 2 C R 3 . The direction of polar- 
ization p lies in the tangent plane TkS 2 . The frame F is fully determined by k 
and p and contains all information concerning the polarization and direction of 
propagation. (The third vector of F is indeed k x p). 

We follow the wave state along a ray using the frame F . Between scatter- 
ing events, F remains constant. The changes of F occur at scattering events. 
As F is a frame, any change, or jump, of F corresponds to a rotation matrix. 
The scattering angle 9 is random, following a probability distribution func- 
tion (pdf) called phase function that depends on a/ A and on the nature 
of the scatterer. The average value of (cos 6) = J <£>(#) cos (9 sin #d# is called 
the scattering anisotropy and is noted g. In the forward scattering regime, the 
main scattering angle 9 is small and g is close to 1. We use the ZYZ conven- 
tion for the Euler angles describing the rotations of SO(3) and note the angles 

(l/j,0,(f) G ]— 7T,7r] X ]0,7r] X ]— 7T,7r]. 

The rotation matrix P acting on the frame F at a scattering event is fully 
determined by the incoming direction k and the outgoing direction k' and the 
requirement that the vector p is parallel transported |15j . On the unit sphere, 
jumps of the stochastic process are represented by geodesies (arcs of great cir- 
cles). Moreover, parallel transport in the direct space implies parallel transport 
in the phase space from TkS 2 to TyS 2 [TT] and therefore that the angle be- 
tween p and the geodesic remains constant. As a consequence the Euler angles 
ip and <p of the rotation matrix P are related by (see figure [TJ 

i, = -<p. (2) 
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Figure 1: Parallel transport of a vector along the geodesic of length 9 (solid 
line). The parallel transport constraint results in the fact that ip and tp are 
equal (algebrically opposite). 



Parallel transport of polarization through a scattering event may therefore 
be described by a rotation matrix P(0,ip), is the length of the geodesic on 
the unit sphere and ip is the angle between this geodesic and the polarization 
vector p. The relation between the incoming and outgoing frames is 



Note that the rotation matrix P acts on the right of F. An expression in 
terms of left action of rotation matrices can be obtained, but would lead to a 
random process over 5*0(3) with dependent increments. The independence of 
increments will be useful for the repeated action |3) on f that we present in 
the next section. 

3 Multiple scattering process 

We consider a plane wave with linear polarization, corresponding to the frame 
Fq as described in the previous section. The scattering events are described as 
random rotations acting on the right of Fq. The frame Fq changes according 
to Equation ((3]) for each scaterring event, so that it becomes the frame F n after n 
scattering events: 



The scattering event are independent and the scatterers are identical, so that 
the random rotations P (0 m , ip m ) are independent and identically distributed. 



F' = FP(9,ip). 



(3) 



F n = F Q P (0 1; Vl) P (02, fo) ■ ■ "P (9n, l/> n ) ■ 



(4) 
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From the source to the observer, a wave may encounter a variable number 
of scatterers. The number of scattering events in a time t can be taken as a 
Poisson process N(t) [T3]. The Poisson parameter r\ is directed related to the 
mean free path £ by -q — c/l. All these considerations result in the process F t 
which we define as 

N(t) 

Ft = F o Y[P(0 k ,Tp k ), (5) 

fe=i 

where denotes the right-sided product on 50(3). This equation expresses 
the stochastic process F t as a compound Poisson process (CPP) on 50(3) with 
parallel transport. We denote by p v the distribution of the process F t at a given 
time t, with v — rjt. A Poisson process is a pure jump process with independent 
increments and we have 

71=0 

where e _l/ ^j is the probability to have exactly n scattering events in the time 
interval [0, i\. Thanks to Equation ((4]) and the independence of the parallel 
transport rotations, pf„ is given by [4]: 

PF n = Pf *PPi*---* PP n = ■ (7) 

Here the symbol * is the convolution product on 50(3) and $*M denotes the 
result of the convolution of n identical functions $. We have used the initial 
condition pp = 5 (I3 — Fq). In the case where the source is a distribution of di- 
rections of propagation and polarization, the distribution p p n is the convolution 
of the initial distribution pp with 



4 Harmonic analysis 

We use harmonic analysis on 50(3) to transform the expression ([7]). Harmonic 
analysis on 50(3) is similar to Fourier series and in particular transforms con- 
volutions into products. As in Ref. |13j . we get a semi-analytic expression of p v 
that we use for statistical estimation in Section [5j 

Th Fourier basis of functions / € L 2 (50(3),R) is made of the Wigner D- 
matrices (D^)j^o 01!]. The Wigner D-matrices are (2j + 1) x (2j + 1) square 
matrices with coefficients 

£>^„(V, 9, <p) = e-^< B (fi)e-^. (8) 

The Fourier coefficients are matrices of the same size which coefficients are de- 
fined as n = (<J> I D 3 m „), where ( ) represents the scalar product for function 
on 50(3) [B]. Like in the circular Fourier theory, the transform of a convolution 
is a product (here, a matrix product) 



We deduce the Fourier coefficients of p„ from Equation ^ 

oo „ 

w — ^ v t~ A n r ~ . i 

S^E^— =exp[K* J -^ + i)J 5 (10) 

n=0 

where exp denotes the matrix exponential. Using the notation A(x) — 2n ^ fc S(x- 
2irk) to put the parallel transport constraint on the Fourier coefficients we ob- 
tain 

$ m,n = ^JJJ HWV + ^D^M, 9, riBmddSdipty. (11) 
We can simplify using the definition ([8j 

*l,n = % j[" $(0)< m (0)sin0d0, (12) 

and show that the Fourier matrices $ 3 are diagonal for all j. Thanks to Equa- 
tion (|10[) . it is clear that p^ 3 is also diagonal at all orders j. 
Using the inverse Fourier formula [BJ, p v reads: 

P,= 2-£(2i+l)TV(^t) (is) 

where Tr denotes the trace and the t denotes the Hermitian conjugation of 
elements from A^2j+i(C). 

The variable f3 = if> — ip represents the direction of polarization of the wave 
after its propagation through the medium, which is the geometric phase. Thanks 
to symmetry, the distribution of /3 is the same as the one of ip + ip. Combined 
with the fact that p^ J is diagonal, we conclude that p v is a function of 9 and j5 

p v (9, 0) = Ro(9, v) + 2 £ cos(m/3)il m (0, v), (14) 

m>l 

R m (9, v) = i- £ (2j + l)e"C5S..m-D4 fm (e). (15) 

Thus, Equation (ffT]) is a general expression of the distribution of geometric 
phase in multiple scattering regime for polarized waves. We have only assumed 
that the scatterers are spherical and identical. 

The limit of Brownian motion, which is made of isotropic infinitesimal steps 
of spherical length 5 occuring at a large rate 77, corresponds to the physical 
situation where scatterers are very weak but have a large spatial density. For 
a single step, we get <& J m m = exp[— hj(J + l)^ 2 ]- Taking the limit S — > 
and 77 — > 00 such that 8 2 v — S 2 r)t — Dt remains constant, formula (fl"4")l gives the 
distribution of the geometric phase for a wave propagating in a heterogeneous 
continuous medium. In formula (fl4l) . one should replace v{& m m — 1) by — 
l)Dt. This formula was demonstrated by Perrin in his study of rotational 
Brownian motion |17[ 118]. while computing the distribution of tp + ip. 
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Henyey-Greenstein 




Figure 2: Evolution of the density {p^ti^ = 0, ft) with time in a medium con- 
taining Henyey-Greenstein scatterers with anisotropy g = 0.8. The curves are 
iso-density levels from 0.05 to 0.45 by increments of 0.05. The vertical lines 
correspond to the slices displayed on Figure [3] 




Figure 3: Distribution of the geometric phase ft for the values of the Poisson 
parameter v = r\t : v = 20, v = 25, v = 30 and v = 35. corresponding to the 
vertical lines of Figure [2j 
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Gaussian 




Dt 



Figure 4: Evolution of the density pDt(& — 0, j3) with time for a Gaussian phase 
function. This corresponds to the rotational Brownian motion on the sphere 
studied by Perrin in Ref. [18]. The curves are iso-density levels from 0.05 to 
0.45 by increments of 0.05. The vertical lines correspond to the slices displayed 
on Figure [5j 




Figure 5: Density pDt(6 = 0, 0) for a Gaussian phase function with the val- 
ues Dt = 0.03, Dt = 0.05, Dt = 0.07 and Dt = 0.09 corresponding to the 
vertical lines of Figure [4] 
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Taking 9 = in Equation (fH|) we get the result demonstrated by Antoine et 
al. [2] , and integrating over 9 we find the result obtained by Krishna et al. [12] up 
to a small difference, has to be replaced by + — m 2 ). Equation (TH|) 

is therefore a generalization of these known formulas that can be applied to an 
arbitrary phase function, also in the case of a small number of scattering events 
and for an arbitrary deviation angle 9. The origin of the small difference comes 
from the fact that the Brownian motion of the tangent vector of a particle mov- 
ing in space R 3 has one more degree of freedom than the Brownian motion on S 2 
(the rotation degree of freedom about the tangent vector). In References [2"1[T2"] 
the extra degree of freedom disappears because these authors considered rota- 
tions with -0 = and in this work, we considered rotations with ip + tp = 0. Let 
us remark that the set of these latter rotations forms a group but that the set 
of rotations with -0 = does not. 

We display the distribution obtained from formula (TH|) in the two following 
cases: discrete Henyey-Greenstein scatterers with anisotropy 5 = 0.8 on a range 
of Poisson parameter up to v — 20 on the one hand (Figures O and [3]) and 
diffusive limit with D — 1 rad 2 s _1 on the other hand (Figures 2] and [5]) . The 
exit angle is 8 = for all of the figures. We can observe that the behavior of 
the distributions on Figures [3] and [5] are very different, especially at short times, 
the distribution p v is sharply peeked for the Henyey-Greenstein scatterers but 
in the case of Brownian motion it is smooth. At long times, depolarization 
is observed in both cases. The depolarized state is reached with significantly 
distinct behaviours. This is an illustration that depolarization is significantly 
dependent on the scattering properties of the medium. It is therefore natural to 
investigate how these scattering properties can be estimated from observations 
through signal analysis. 

We have shown that the polarization distribution depends on the angle be- 
tween the outgoing wave vector fe out and the initial one fco. Moreover, it ac- 
tually depends on the position of the last scattering event. Let us call \ the 
angle coordinate of this point in the cylindrical coordinates of axis fco in the 
frame F , and (9, 4>) the spherical coordinates of fc out in the same frame. The 
distribution of polarization is then the one given by formula (| 14|) shifted by an 
angle (1 — cos 9)x- This shift is due to the fact that the reference j3 = is taken 
in the parallel transported frame from the source to the last scattering event, so 
that an extra phase appears in the frame F b ■ This extra phase has been already 
observed in backscattering configuration (9 = ir) [8j |2"T]. 

We remark that in the case where m depends only on j (and is noted <fr J ), 
formula (|14[) simplifies to 

M9, P) = J> + i^-D ^ + lh , (16) 
rrf. sin 2^ 

with cos ^ = cos I cos ^. As a consequence, there is an unexpected symmetry 
for p v between the outgoing angle and the geometric phase. 
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5 Inverse problem 

In this section, we take advantage of the Fourier expansion and propose a statis- 
tical estimation of the mean free path I obtained from the measured distribution 
of the geometric phase j3. This is achieved by estimating the Poisson parame- 
ter v = rjt. The estimate of v is denoted by v . 

The inverse problem is solved using an expectation-minimization (EM) ap- 
proach 5 with the parametric CPP. The EM algorithm is based on the max- 
imization of the log-likelihood of the a posteriori distribution given in Equa- 
tion ((6|) . Suppose that we are given a sample of size M of observations (mea- 
surements) F m , 1 ^ to ^ M. We assume that each observation F m follows 
independently the probability law (fl4|) . the joint probability distribution func- 
tion is simply the product Y\ m Pu{Fm), with p v {F) = p v (9,ip + ip), where ip, 
9 and ip are the Euler angles of F . The value of v that maximizes this log- 
likelihood is 



The EM algorithm is an iterative procedure that almost surely converges 
to a local minimum [5 . The minimization step (M-step) consists in updat- 
ing the estimate C>, while the expectation step (E-step) consists in updating 
P (n | F, Vi) — pp n (F)/pc, i (F) x e~ Vi vf/n\. We get a sequence of estimates 
with the relation 



In practice the sums over n in Equation (|18l) and Equation ([5]) and over j in 
Equation (TT5)) are truncated to an arbitrarily fixed value N. At most N scat- 
tering events are described by the model, therefore N should be taken large 
compared to L/£. This EM-algorithm is able to estimate one parameter, the 
Poisson parameter v (or equivalently the mean free path I — ct/v), which means 
that we have considered that the phase function $ is known. It is nonetheless 
possible to estimate simultaneously any Unite number of parameters of p„, such 
as the anisotropy g of the scatterers, but the iteration formulas are more so- 
phisticated. 

We illustrate the behaviour of the EM algorithm in Figure [6] which shows 
the convergence of the estimator v for different initial values vq. Acceleration 
procedures could be developped in the case of large dataset from experimental 
setup. Note that no local minima were reached in the simulation, leading to 
convergence in all the cases. However, as can be seen for initial values inferior 
to 10 in Figure El underestimated initial values lead to a systematic bias in v. 
This suggests that high values should be privileged when processing dataset, as 
it does not penalize convergence and leads to a more accurate estimate. 




(17) 




(18) 
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Figure 6: Convergence of the estimate of the Poisson parameter rji using EM 
algorithm with different initialization values ranging from 1 to 29. The actual 
value of rj was set to 10. 



6 Conclusions 

The description of the depolarization of multiply scattered waves can be made, 
in the forward scattering regime, through a compound Poisson process (CPP). 
This stochastic model predicts the distribution of the geometric phase in all 
directions, for discrete or continuous scattering media, generalizing existing re- 
sults (forward outgoing direction or a spatially fluctuating medium). Moreover, 
the CPP model allows a more detailed description of the phenomenon as it 
provides the behaviour of this distribution as a function of the output scatter- 
ing angle 9. The present approach allows to design an iterative procedure to 
estimate properties of the scattering medium through the measurement of the 
outgoing polarization distribution. We have illustrated this point by presenting 
an expectation-minimization (EM) algorithm for the estimation of the Poisson 
parameter, which is directly linked to the transport mean free path £*. An in- 
teresting feature of this technique is that is relies on polarization rather than on 
amplitude measurements. Future work will consist in validating the proposed 
model and estimation algorithm on measurement of polarization distribution 
for real experiment. 
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